1/f Noise Project Manuscript 1 Final Analysis
Fine-Scale Spatial Heterogeneity in Phytochemistry Increases Herbivore Performance by Enhancing Herbivore Foraging
Vincent S. Pan, Enakshi Ghosh, Paul J. Ode, William C. Wetzel, Kadeem J. Gilbert, Ian S. Pearse
May 17, 2024
1 Setup
1.1 Load libraries
library(survival)
library(piecewiseSEM)
source("spat1f/init_analysis.R")Data last updated on 2024-05-031.2 Handy functions not in spat1f
trend_by_pre_wt <- function(model, term){
emmeans::emtrends(model, var = "cat_pre_wt_log_scale", specs = term) %>%
pairs()
}
contrast_by_pre_wt <- function(model, term, type = c("pairwise", "trt.vs.ctrl"), at = c(-2, 2)){
type <- match.arg(type)
f <- as.formula(sprintf("%s ~ %s | cat_pre_wt_log_scale", type, term))
emmeans::emmeans(model,
f,
var = term,
at = list(cat_pre_wt_log_scale = at))
}
check_mod <- function(model){
DHARMa::simulateResiduals(model) %>%
plot()
performance::posterior_predictive_check(model)
}1.3 Quick cleaning
# Main data.frame for analyses
d <- ref_data %>%
filter(error == 0) %>%
mutate(
cat_pre_wt_log = log(cat_pre_wt),
cat_pre_wt_log_scale = as.numeric(scale(log(cat_pre_wt))),
cat_size = ifelse(cat_pre_wt < median(cat_pre_wt), "small", "big"),
has_var = ifelse(var_trt != "constant", 1, 0),
mean_toxic_conc_scale = as.numeric(scale(mean_toxic_conc)),
mean_toxic_scale = as.numeric(scale(mean_toxic)),
area_herb_log_scale = as.numeric(scale(log(area_herb+1))),
var_toxic_12_scale = as.numeric(scale(var_toxic_12)),
on_toxic_logit_scale = as.numeric(scale(adjust_prop(on_toxic,
trans = "emp",
nudge.method = "none",
na.action = "ignore"))),
select_scale = as.numeric(scale(less_toxic)),
ava_qual_scale = as.numeric(scale(ava_qual)),
ava_qual_logit_scale = as.numeric(scale(adjust_prop(ava_qual,
trans = "emp",
nudge.method = "none",
na.action = "ignore"))),
sl_mean_obs_log_scale = as.numeric(scale(log(sl_mean_obs))),
cat_pre_wt_log_scale_sq = as.numeric(scale(log(cat_pre_wt)^2)),
RGR_scale = as.numeric(scale(RGR)),
sl_mean_pred = shape * scale
)
# Subset of data.frame for SEM
d2 <- d %>%
filter(
var_trt != "constant"
) %>%
filter(
!is.na(area_herb_log_scale) &
!is.na(var_toxic_12_scale) &
!is.na(mean_toxic_conc_scale) &
!is.na(on_toxic_logit_scale) &
!is.na(select_scale) &
!is.na(ava_qual_scale) &
!is.na(sl_mean_obs) &
!is.na(RGR)) %>%
mutate(
var_high = ifelse(var_trt == "high_var", 1, 0),
beta_red = ifelse(as.character(beta) == 5, 1, 0),
beta_white = ifelse(as.character(beta) == 0, 1, 0),
beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
)2 Herbivore Performance
2.1 RGR
2.1.1 Model summary
rgr_m <- glmmTMB(
RGR ~
(var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # residual plot show significant non-linearity
(1|session_id),
data = d %>%
filter(
var_trt != "constant"
)
); summary(rgr_m) Family: gaussian ( identity )
Formula:
RGR ~ (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
(1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
-879.2 -849.4 450.6 -901.2 100
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.979e-06 0.001995
Residual 1.603e-05 0.004004
Number of obs: 111, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.6e-05
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0162504 0.0013457 12.075 < 2e-16 ***
var_trtlow_var -0.0034872 0.0007740 -4.506 6.62e-06 ***
beta0 -0.0003065 0.0009664 -0.317 0.75112
beta5 -0.0013545 0.0009173 -1.477 0.13978
cat_pre_wt_log_scale -0.0048759 0.0011991 -4.066 4.78e-05 ***
I(cat_pre_wt_log_scale^2) -0.0021953 0.0006731 -3.262 0.00111 **
var_trtlow_var:cat_pre_wt_log_scale 0.0025873 0.0007897 3.276 0.00105 **
beta0:cat_pre_wt_log_scale 0.0011924 0.0009491 1.256 0.20897
beta5:cat_pre_wt_log_scale 0.0024731 0.0009710 2.547 0.01087 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(rgr_m, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: RGR
Chisq Df Pr(>Chisq)
(Intercept) 145.8141 1 < 2.2e-16 ***
var_trt 20.3008 1 6.617e-06 ***
beta 2.3963 2 0.301758
cat_pre_wt_log_scale 16.5348 1 4.777e-05 ***
I(cat_pre_wt_log_scale^2) 10.6374 1 0.001108 **
var_trt:cat_pre_wt_log_scale 10.7334 1 0.001052 **
beta:cat_pre_wt_log_scale 6.4936 2 0.038899 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.1.2 Post-hoc contrasts
trend_by_pre_wt(rgr_m, "var_trt") contrast estimate SE df t.ratio p.value
high_var - low_var -0.00259 0.00079 100 -3.276 0.0014
Results are averaged over the levels of: beta contrast_by_pre_wt(rgr_m, "var_trt")$emmeans
cat_pre_wt_log_scale = -2:
var_trt emmean SE df lower.CL upper.CL
high_var 0.014224 0.00355 100 0.007178 0.02127
low_var 0.005562 0.00312 100 -0.000635 0.01176
cat_pre_wt_log_scale = 2:
var_trt emmean SE df lower.CL upper.CL
high_var -0.000393 0.00208 100 -0.004512 0.00373
low_var 0.001295 0.00238 100 -0.003421 0.00601
Results are averaged over the levels of: beta
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
high_var - low_var 0.00866 0.00181 100 4.790 <.0001
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
high_var - low_var -0.00169 0.00171 100 -0.988 0.3256
Results are averaged over the levels of: beta trend_by_pre_wt(rgr_m, "beta") contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.00119 0.000949 100 -1.256 0.4232
(beta-5) - beta5 -0.00247 0.000971 100 -2.547 0.0329
beta0 - beta5 -0.00128 0.000940 100 -1.363 0.3644
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates contrast_by_pre_wt(rgr_m, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df lower.CL upper.CL
-5 0.012890 0.00382 100 0.005313 0.02047
0 0.010199 0.00333 100 0.003591 0.01681
5 0.006589 0.00319 100 0.000255 0.01292
cat_pre_wt_log_scale = 2:
beta emmean SE df lower.CL upper.CL
-5 -0.001439 0.00226 100 -0.005922 0.00304
0 0.000639 0.00231 100 -0.003940 0.00522
5 0.002153 0.00257 100 -0.002946 0.00725
Results are averaged over the levels of: var_trt
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 0.00269 0.00221 100 1.218 0.4456
(beta-5) - beta5 0.00630 0.00223 100 2.820 0.0158
beta0 - beta5 0.00361 0.00211 100 1.711 0.2062
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.00208 0.00205 100 -1.015 0.5689
(beta-5) - beta5 -0.00359 0.00206 100 -1.746 0.1935
beta0 - beta5 -0.00151 0.00211 100 -0.719 0.7529
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates 2.1.3 Misc stats
r2(rgr_m, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.555
Marginal R2: 0.4442.1.4 Check model
check_mod(rgr_m)
2.1.5 Plot
g1 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) +
geom_line(aes(color = beta), linewidth = 2) +
geom_point(
data = rgr_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = beta),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill","color")) +
labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
color = expression(beta), fill = expression(beta))
g2 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","var_trt")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = var_trt), alpha = 0.2) +
geom_line(aes(color = var_trt), linewidth = 2) +
geom_point(
data = rgr_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = var_trt),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
theme(legend.position = "top") +
scale_color_discrete(type = c("navy","skyblue"),
aesthetics = c("fill","color"),
label = c("High", "Low")) +
labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
color = "Variation", fill = "Variation")
ggarrange(g1, g2 + labs(y = ""))
2.2 Time to pupation
2.2.1 Model summary
pupt_m_full <- glmmTMB(
time_to_pupation ~
(var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = Gamma(link = "log"),
); summary(pupt_m_full) Family: Gamma ( log )
Formula: time_to_pupation ~ (var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
338.1 366.3 -158.0 316.1 85
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.839e-11 5.328e-06
Number of obs: 96, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0254
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.180302 0.036853 59.16 < 2e-16 ***
var_trtlow_var 0.102145 0.033998 3.00 0.00266 **
beta0 -0.001845 0.040823 -0.05 0.96395
beta5 0.059543 0.041181 1.45 0.14821
cat_pre_wt_log_scale -0.428662 0.032643 -13.13 < 2e-16 ***
I(cat_pre_wt_log_scale^2) -0.081267 0.019069 -4.26 2.03e-05 ***
var_trtlow_var:cat_pre_wt_log_scale 0.026518 0.038876 0.68 0.49516
beta0:cat_pre_wt_log_scale -0.030363 0.043195 -0.70 0.48210
beta5:cat_pre_wt_log_scale -0.126801 0.047095 -2.69 0.00709 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(pupt_m_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: time_to_pupation
Chisq Df Pr(>Chisq)
(Intercept) 3500.2207 1 < 2.2e-16 ***
var_trt 9.0264 1 0.002661 **
beta 2.7351 2 0.254736
cat_pre_wt_log_scale 172.4448 1 < 2.2e-16 ***
I(cat_pre_wt_log_scale^2) 18.1614 1 2.029e-05 ***
var_trt:cat_pre_wt_log_scale 0.4653 1 0.495162
beta:cat_pre_wt_log_scale 7.5589 2 0.022836 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pupt_m <- glmmTMB(
time_to_pupation ~
var_trt + beta * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = Gamma(link = "log"),
); summary(pupt_m) Family: Gamma ( log )
Formula:
time_to_pupation ~ var_trt + beta * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
(1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
336.5 362.2 -158.3 316.5 86
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.159e-11 4.647e-06
Number of obs: 96, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0255
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.180555 0.036937 59.03 < 2e-16 ***
var_trtlow_var 0.107447 0.033197 3.24 0.00121 **
beta0 -0.003578 0.040842 -0.09 0.93019
beta5 0.059536 0.041271 1.44 0.14915
cat_pre_wt_log_scale -0.420147 0.030221 -13.90 < 2e-16 ***
I(cat_pre_wt_log_scale^2) -0.083365 0.018881 -4.42 1.01e-05 ***
beta0:cat_pre_wt_log_scale -0.026011 0.042849 -0.61 0.54383
beta5:cat_pre_wt_log_scale -0.126834 0.047183 -2.69 0.00719 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(pupt_m, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: time_to_pupation
Chisq Df Pr(>Chisq)
(Intercept) 3485.0791 1 < 2.2e-16 ***
var_trt 10.4759 1 0.001209 **
beta 2.8053 2 0.245947
cat_pre_wt_log_scale 193.2743 1 < 2.2e-16 ***
I(cat_pre_wt_log_scale^2) 19.4946 1 1.009e-05 ***
beta:cat_pre_wt_log_scale 7.6890 2 0.021397 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2.2 Post-hoc contrasts
contrast_by_pre_wt(pupt_m, "var_trt")$emmeans
cat_pre_wt_log_scale = -2:
var_trt emmean SE df asymp.LCL asymp.UCL
high_var 2.808 0.0769 Inf 2.657 2.96
low_var 2.915 0.0818 Inf 2.755 3.08
cat_pre_wt_log_scale = 2:
var_trt emmean SE df asymp.LCL asymp.UCL
high_var 0.924 0.0672 Inf 0.792 1.06
low_var 1.031 0.0725 Inf 0.889 1.17
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df z.ratio p.value
high_var - low_var -0.107 0.0332 Inf -3.237 0.0012
cat_pre_wt_log_scale = 2:
contrast estimate SE df z.ratio p.value
high_var - low_var -0.107 0.0332 Inf -3.237 0.0012
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale. trend_by_pre_wt(pupt_m, "beta") contrast estimate SE df z.ratio p.value
(beta-5) - beta0 0.026 0.0428 Inf 0.607 0.8163
(beta-5) - beta5 0.127 0.0472 Inf 2.688 0.0197
beta0 - beta5 0.101 0.0473 Inf 2.132 0.0834
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates contrast_by_pre_wt(pupt_m, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df asymp.LCL asymp.UCL
-5 2.741 0.0941 Inf 2.557 2.93
0 2.790 0.0971 Inf 2.599 2.98
5 3.054 0.1088 Inf 2.841 3.27
cat_pre_wt_log_scale = 2:
beta emmean SE df asymp.LCL asymp.UCL
-5 1.061 0.0845 Inf 0.895 1.23
0 1.005 0.0805 Inf 0.847 1.16
5 0.866 0.0917 Inf 0.687 1.05
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df z.ratio p.value
(beta-5) - beta0 -0.0484 0.1018 Inf -0.476 0.8828
(beta-5) - beta5 -0.3132 0.1133 Inf -2.763 0.0158
beta0 - beta5 -0.2648 0.1136 Inf -2.331 0.0516
cat_pre_wt_log_scale = 2:
contrast estimate SE df z.ratio p.value
(beta-5) - beta0 0.0556 0.0875 Inf 0.636 0.8005
(beta-5) - beta5 0.1941 0.0915 Inf 2.122 0.0854
beta0 - beta5 0.1385 0.0928 Inf 1.492 0.2945
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates 2.2.3 Misc stats
r2(pupt_m, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.8812.2.4 Check model
check_mod(pupt_m)
2.2.5 Plot
g3 <- marginal_effects(pupt_m, terms = c("cat_pre_wt_log_scale","beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) +
geom_line(aes(color = beta), linewidth = 2) +
geom_point(
data = pupt_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = time_to_pupation, color = beta),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log2") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill","color")) +
labs(x = "Cat pre-weight (g)", y = "Time to pupation (days)",
color = expression(beta), fill = expression(beta))
g4 <- marginal_effects(pupt_m, terms = c("var_trt")) %>%
ggplot(aes(x = var_trt, y = yhat)) +
geom_point(
data = pupt_m$frame,
aes(x = var_trt, y = time_to_pupation, color = var_trt),
position = position_jitter(width = 0.2),
size = 3, alpha = 0.8
) +
geom_pointrange(
aes(ymax = upper, ymin = lower),
color = "black",
size = 1, linewidth = 2, shape = 1,
) +
theme_bw(base_size = 15) +
scale_y_continuous(trans = "log2") +
theme(legend.position = "top") +
scale_color_discrete(type = c("navy","skyblue"),
aesthetics = c("color"),
label = c("High", "Low")) +
labs(x = "Variation", y = "Time to pupation (days)",
color = "Variation")
ggarrange(g3, g4 + labs(y = ""))
2.3 Eclosure
2.3.1 Model summary
eclosure_m_full <- glmmTMB(
eclosed ~
(var_trt + beta) * cat_pre_wt_log_scale +
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = binomial(link = "logit"),
); summary(eclosure_m_full) Family: binomial ( logit )
Formula:
eclosed ~ (var_trt + beta) * cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
185.3 210.9 -83.6 167.3 119
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.1661 0.4076
Number of obs: 128, groups: session_id, 5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.53530 0.41050 1.304 0.192
var_trtlow_var 0.01707 0.37507 0.046 0.964
beta0 -0.21544 0.46265 -0.466 0.641
beta5 -0.50749 0.45596 -1.113 0.266
cat_pre_wt_log_scale 0.10452 0.36089 0.290 0.772
var_trtlow_var:cat_pre_wt_log_scale 0.17526 0.36489 0.480 0.631
beta0:cat_pre_wt_log_scale 0.24778 0.43553 0.569 0.569
beta5:cat_pre_wt_log_scale 0.20159 0.44699 0.451 0.652car::Anova(eclosure_m_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: eclosed
Chisq Df Pr(>Chisq)
(Intercept) 1.7004 1 0.1922
var_trt 0.0021 1 0.9637
beta 1.2502 2 0.5352
cat_pre_wt_log_scale 0.0839 1 0.7721
var_trt:cat_pre_wt_log_scale 0.2307 1 0.6310
beta:cat_pre_wt_log_scale 0.3659 2 0.8328eclosure_m <- glmmTMB(
eclosed ~
(var_trt + beta) + cat_pre_wt_log_scale +
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = binomial(link = "logit"),
); summary(eclosure_m) Family: binomial ( logit )
Formula:
eclosed ~ (var_trt + beta) + cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
179.9 197.0 -84.0 167.9 122
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.1676 0.4094
Number of obs: 128, groups: session_id, 5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.551286 0.413578 1.333 0.183
var_trtlow_var 0.005064 0.373402 0.014 0.989
beta0 -0.237692 0.460689 -0.516 0.606
beta5 -0.523521 0.456587 -1.147 0.252
cat_pre_wt_log_scale 0.330169 0.220049 1.500 0.134car::Anova(eclosure_m, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: eclosed
Chisq Df Pr(>Chisq)
var_trt 0.0002 1 0.9892
beta 1.3212 2 0.5165
cat_pre_wt_log_scale 2.2513 1 0.13352.3.2 Post-hoc contrasts
2.3.3 Misc stats
r2(eclosure_m, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.092
Marginal R2: 0.0462.3.4 Check model
check_mod(eclosure_m)
2.3.5 Plot
2.4 Survival
2.4.1 Model summary
surv_m_full <- coxph(
Surv(surv_time, observed_dead) ~
(var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
strata(session_id),
data = d %>%
filter(var_trt != "constant") %>%
mutate(
beta = as.factor(as.character(beta))
),
); summary(surv_m_full)Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) *
cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id),
data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))
n= 128, number of events= 105
coef exp(coef) se(coef) z Pr(>|z|)
var_trtlow_var -0.1115 0.8945 0.2124 -0.525 0.599727
beta0 -0.0896 0.9143 0.2579 -0.347 0.728306
beta5 -0.2521 0.7772 0.2605 -0.968 0.333144
cat_pre_wt_log_scale 0.9641 2.6224 0.2556 3.772 0.000162
I(cat_pre_wt_log_scale^2) 0.2450 1.2776 0.1447 1.693 0.090408
var_trtlow_var:cat_pre_wt_log_scale -0.1626 0.8499 0.1998 -0.814 0.415661
beta0:cat_pre_wt_log_scale -0.2230 0.8001 0.2309 -0.966 0.334053
beta5:cat_pre_wt_log_scale -0.4789 0.6195 0.2584 -1.854 0.063807
var_trtlow_var
beta0
beta5
cat_pre_wt_log_scale ***
I(cat_pre_wt_log_scale^2) .
var_trtlow_var:cat_pre_wt_log_scale
beta0:cat_pre_wt_log_scale
beta5:cat_pre_wt_log_scale .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var 0.8945 1.1179 0.5899 1.356
beta0 0.9143 1.0937 0.5515 1.516
beta5 0.7772 1.2867 0.4665 1.295
cat_pre_wt_log_scale 2.6224 0.3813 1.5892 4.328
I(cat_pre_wt_log_scale^2) 1.2776 0.7827 0.9621 1.696
var_trtlow_var:cat_pre_wt_log_scale 0.8499 1.1766 0.5746 1.257
beta0:cat_pre_wt_log_scale 0.8001 1.2499 0.5089 1.258
beta5:cat_pre_wt_log_scale 0.6195 1.6142 0.3734 1.028
Concordance= 0.637 (se = 0.047 )
Likelihood ratio test= 22 on 8 df, p=0.005
Wald test = 22.31 on 8 df, p=0.004
Score (logrank) test = 24.57 on 8 df, p=0.002drop1(surv_m_full, test = "Chisq")Single term deletions
Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + strata(session_id)
Df AIC LRT Pr(>Chi)
<none> 459.50
I(cat_pre_wt_log_scale^2) 1 460.40 2.90 0.08841 .
strata(session_id) 0 787.87 328.37
var_trt:cat_pre_wt_log_scale 1 458.16 0.67 0.41395
beta:cat_pre_wt_log_scale 2 458.98 3.48 0.17532
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1surv_m <- coxph(
Surv(surv_time, observed_dead) ~
(var_trt + beta) + cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
strata(session_id),
data = d %>%
filter(var_trt != "constant") %>%
mutate(
beta = as.factor(as.character(beta))
),
); summary(surv_m)Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) +
cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id),
data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))
n= 128, number of events= 105
coef exp(coef) se(coef) z Pr(>|z|)
var_trtlow_var -0.07077 0.93167 0.20845 -0.340 0.734212
beta0 -0.07410 0.92858 0.25808 -0.287 0.774026
beta5 -0.17722 0.83759 0.25804 -0.687 0.492206
cat_pre_wt_log_scale 0.65947 1.93376 0.18164 3.631 0.000283 ***
I(cat_pre_wt_log_scale^2) 0.29584 1.34425 0.14124 2.095 0.036205 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var 0.9317 1.0733 0.6192 1.402
beta0 0.9286 1.0769 0.5599 1.540
beta5 0.8376 1.1939 0.5051 1.389
cat_pre_wt_log_scale 1.9338 0.5171 1.3545 2.761
I(cat_pre_wt_log_scale^2) 1.3443 0.7439 1.0192 1.773
Concordance= 0.602 (se = 0.047 )
Likelihood ratio test= 18.08 on 5 df, p=0.003
Wald test = 18.29 on 5 df, p=0.003
Score (logrank) test = 19.42 on 5 df, p=0.002drop1(surv_m, test = "Chisq")Single term deletions
Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) + cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + strata(session_id)
Df AIC LRT Pr(>Chi)
<none> 457.42
var_trt 1 455.53 0.12 0.7341150
beta 2 453.89 0.48 0.7883535
cat_pre_wt_log_scale 1 469.05 13.63 0.0002226 ***
I(cat_pre_wt_log_scale^2) 1 459.84 4.42 0.0355376 *
strata(session_id) 0 784.70 327.28
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.4.2 Post-hoc contrasts
2.4.3 Misc stats
r2(surv_m, tolerance = 10^-10) Nagelkerke's R2: 0.1352.4.4 Check model
cox.zph(surv_m) chisq df p
var_trt 1.60 1 0.21
beta 1.43 2 0.49
cat_pre_wt_log_scale 2.61 1 0.11
I(cat_pre_wt_log_scale^2) 1.05 1 0.31
GLOBAL 5.53 5 0.352.4.5 Plot
3 Structural Equation Model
3.1 Sub-models
3.1.1 RGR
subm_rgr <- glmmTMB(
RGR_scale ~
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
mean_toxic_conc_scale +
area_herb_log_scale * cat_pre_wt_log_scale +
sl_mean_obs_log_scale +
var_toxic_12_scale +
(1|session_id),
data = d2
); summary(subm_rgr) Family: gaussian ( identity )
Formula: RGR_scale ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +
mean_toxic_conc_scale + area_herb_log_scale * cat_pre_wt_log_scale +
sl_mean_obs_log_scale + var_toxic_12_scale + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
140.2 165.6 -60.1 120.2 84
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01202 0.1096
Residual 0.20224 0.4497
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.202
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.041726 0.078479 0.532 0.59494
cat_pre_wt_log_scale -3.367086 0.405069 -8.312 < 2e-16
cat_pre_wt_log_scale_sq -2.789806 0.417379 -6.684 2.32e-11
mean_toxic_conc_scale -0.233908 0.068928 -3.393 0.00069
area_herb_log_scale 0.773436 0.091519 8.451 < 2e-16
sl_mean_obs_log_scale -0.164921 0.053712 -3.070 0.00214
var_toxic_12_scale -0.005495 0.055009 -0.100 0.92043
cat_pre_wt_log_scale:area_herb_log_scale -0.200719 0.091201 -2.201 0.02775
(Intercept)
cat_pre_wt_log_scale ***
cat_pre_wt_log_scale_sq ***
mean_toxic_conc_scale ***
area_herb_log_scale ***
sl_mean_obs_log_scale **
var_toxic_12_scale
cat_pre_wt_log_scale:area_herb_log_scale *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.2 Variance in toxin (12 hour)
subm_var_toxic <- glmmTMB(
var_toxic_12_scale ~
beta_numeric_scale +
var_high +
sl_mean_obs_log_scale +
select_scale +
ava_qual_logit_scale +
(1|session_id),
data = d2
); summary(subm_var_toxic) Family: gaussian ( identity )
Formula:
var_toxic_12_scale ~ beta_numeric_scale + var_high + sl_mean_obs_log_scale +
select_scale + ava_qual_logit_scale + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
163.9 184.2 -73.9 147.9 86
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01202 0.1097
Residual 0.27355 0.5230
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.274
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.86243 0.09523 -9.057 < 2e-16 ***
beta_numeric_scale -0.10881 0.05830 -1.866 0.062004 .
var_high 1.69712 0.11500 14.757 < 2e-16 ***
sl_mean_obs_log_scale 0.13145 0.06474 2.030 0.042331 *
select_scale -0.20633 0.08159 -2.529 0.011440 *
ava_qual_logit_scale -0.28324 0.07880 -3.594 0.000325 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.3 Mean step length
subm_sl <- glmmTMB(
sl_mean_obs_log_scale ~
beta_numeric_scale +
on_toxic_logit_scale +
var_high +
select_scale +
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
(1|session_id),
data = d2,
); summary(subm_sl) Family: gaussian ( identity )
Formula:
sl_mean_obs_log_scale ~ beta_numeric_scale + on_toxic_logit_scale +
var_high + select_scale + cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
245.0 267.9 -113.5 227.0 85
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.828e-10 1.682e-05
Residual 6.551e-01 8.094e-01
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.655
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.24784 0.12625 1.963 0.049627 *
beta_numeric_scale 0.18760 0.08691 2.159 0.030879 *
on_toxic_logit_scale 0.38465 0.11614 3.312 0.000926 ***
var_high -0.24283 0.17679 -1.374 0.169587
select_scale 0.35697 0.13067 2.732 0.006297 **
cat_pre_wt_log_scale 1.51783 0.58084 2.613 0.008971 **
cat_pre_wt_log_scale_sq 1.84748 0.61197 3.019 0.002537 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.4 Mean toxin ingested
subm_toxin_ingested <- glmmTMB(
mean_toxic_conc_scale ~
beta_numeric_scale +
var_high * cat_pre_wt_log_scale +
on_toxic_logit_scale +
cat_pre_wt_log_scale_sq +
(1|session_id),
data = d2
); summary(subm_toxin_ingested) Family: gaussian ( identity )
Formula:
mean_toxic_conc_scale ~ beta_numeric_scale + var_high * cat_pre_wt_log_scale +
on_toxic_logit_scale + cat_pre_wt_log_scale_sq + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
154.8 177.7 -68.4 136.8 85
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.90e-11 8.307e-06
Residual 2.51e-01 5.010e-01
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.251
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.26244 0.07904 3.320 0.000899 ***
beta_numeric_scale -0.08727 0.05331 -1.637 0.101636
var_high -0.81255 0.11198 -7.257 3.97e-13 ***
cat_pre_wt_log_scale 0.96560 0.37021 2.608 0.009100 **
on_toxic_logit_scale 0.23369 0.06094 3.835 0.000126 ***
cat_pre_wt_log_scale_sq 1.06006 0.37862 2.800 0.005113 **
var_high:cat_pre_wt_log_scale -0.26558 0.11647 -2.280 0.022596 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.5 Mean time on toxin
subm_on_toxic <- glmmTMB(
on_toxic_logit_scale ~
select_scale + ava_qual_logit_scale +
(1|session_id),
family = gaussian(),
data = d2
); summary(subm_on_toxic) Family: gaussian ( identity )
Formula: on_toxic_logit_scale ~ select_scale + ava_qual_logit_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
-106.0 -93.3 58.0 -116.0 89
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 9.663e-12 3.109e-06
Residual 1.705e-02 1.306e-01
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.0171
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.009281 0.013505 -0.69 0.492
select_scale -0.157732 0.019369 -8.14 3.84e-16 ***
ava_qual_logit_scale -0.909758 0.016169 -56.27 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.6 Area consumption
subm_herb <- glmmTMB(
area_herb_log_scale ~
beta_numeric_scale + var_high * cat_pre_wt_log_scale +
(1|session_id),
data = d2,
); summary(subm_herb) Family: gaussian ( identity )
Formula:
area_herb_log_scale ~ beta_numeric_scale + var_high * cat_pre_wt_log_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
146.4 164.2 -66.2 132.4 87
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.742e-11 8.211e-06
Residual 2.394e-01 4.893e-01
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.239
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.046979 0.074660 0.629 0.52919
beta_numeric_scale 0.001622 0.050860 0.032 0.97457
var_high 0.283395 0.104609 2.709 0.00675 **
cat_pre_wt_log_scale 0.550474 0.087850 6.266 3.7e-10 ***
var_high:cat_pre_wt_log_scale -0.355786 0.113217 -3.143 0.00168 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.7 Toxin selection
subm_select <- glmmTMB(
select_scale ~
beta_numeric_scale + var_high * cat_pre_wt_log_scale +
ava_qual_logit_scale +
(1|session_id),
data = d2,
); summary(subm_select) Family: gaussian ( identity )
Formula:
select_scale ~ beta_numeric_scale + var_high * cat_pre_wt_log_scale +
ava_qual_logit_scale + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
205.9 226.3 -95.0 189.9 86
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 1.243e-10 1.115e-05
Residual 4.415e-01 6.645e-01
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.442
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03222 0.10287 -0.313 0.7541
beta_numeric_scale -0.08040 0.07104 -1.132 0.2577
var_high -0.01374 0.14649 -0.094 0.9253
cat_pre_wt_log_scale -0.12953 0.12030 -1.077 0.2816
ava_qual_logit_scale 0.37617 0.07942 4.736 2.18e-06 ***
var_high:cat_pre_wt_log_scale 0.37178 0.15379 2.417 0.0156 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.8 Neighborhood diet quality
subm_ava <- glmmTMB(
ava_qual_logit_scale ~
var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1|session_id),
family = gaussian(),
data = d2,
); summary(subm_ava) Family: gaussian ( identity )
Formula:
ava_qual_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
248.3 266.1 -117.1 234.3 87
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01866 0.1366
Residual 0.69278 0.8323
Number of obs: 94, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.693
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.20763 0.14063 -1.476 0.13983
var_high 0.45433 0.17311 2.625 0.00868 **
beta_numeric_scale 0.15703 0.09021 1.741 0.08173 .
cat_pre_wt_log_scale 0.18432 0.12145 1.518 0.12912
beta_numeric_scale:cat_pre_wt_log_scale 0.21360 0.09834 2.172 0.02985 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 d-separation test
sem_fit <- psem(
subm_rgr,
subm_var_toxic,
subm_on_toxic,
subm_toxin_ingested,
subm_select,
subm_ava,
subm_herb,
subm_sl,
var_toxic_12_scale %~~% on_toxic_logit_scale, # correlated errors
data = d2
)
sem_summary <- suppressWarnings(suppressMessages(
summary_psem(sem_fit,
no_standardize_x = c("var_high","beta_red","beta_white"),
.progressBar = FALSE,
# Remove 'ava_qual_logit_scale' from the independence claim when 'on_toxic_logit_scale' is conditioned. 'on_toxic_logit_scale' is a more proximal predictor and is highly collinear with ava_qual_logit_scale.
basis.set = keep(basisSet(sem_fit), function(x){
res1 <- (basisSet_predictor(x) == "ava_qual_logit_scale")
res2 <- ("on_toxic_logit_scale" %in% basisSet_conditioned(x))
res <- res1 & res2
return(!res)
}))
))Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.sem_summary$Cstat Fisher.C df P.Value
1 48.845 48 0.439knitr::kable(sem_summary$dTable)| Independ.Claim | Test.Type | DF | Crit.Value | P.Value | |
|---|---|---|---|---|---|
| var_toxic_12_scale ~ cat_pre_wt_log_scale + … | coef | 94 | 1.5236 | 0.1276 | |
| on_toxic_logit_scale ~ cat_pre_wt_log_scale + … | coef | 94 | 0.1994 | 0.8419 | |
| area_herb_log_scale ~ cat_pre_wt_log_scale_sq + … | coef | 94 | 0.8087 | 0.4187 | |
| var_toxic_12_scale ~ cat_pre_wt_log_scale_sq + … | coef | 94 | -1.2587 | 0.2081 | |
| select_scale ~ cat_pre_wt_log_scale_sq + … | coef | 94 | -0.8184 | 0.4131 | |
| ava_qual_logit_scale ~ cat_pre_wt_log_scale_sq + … | coef | 94 | -0.4761 | 0.6340 | |
| on_toxic_logit_scale ~ cat_pre_wt_log_scale_sq + … | coef | 94 | -0.3624 | 0.7171 | |
| RGR_scale ~ beta_numeric_scale + … | coef | 94 | -1.8383 | 0.0660 | |
| on_toxic_logit_scale ~ beta_numeric_scale + … | coef | 94 | 0.1477 | 0.8826 | |
| RGR_scale ~ var_high + … | coef | 94 | -0.3104 | 0.7563 | |
| on_toxic_logit_scale ~ var_high + … | coef | 94 | -0.6569 | 0.5113 | |
| select_scale ~ RGR_scale + … | coef | 94 | -1.0670 | 0.2860 | |
| ava_qual_logit_scale ~ RGR_scale + … | coef | 94 | 1.6016 | 0.1093 | |
| on_toxic_logit_scale ~ RGR_scale + … | coef | 94 | -0.7917 | 0.4285 | |
| sl_mean_obs_log_scale ~ area_herb_log_scale + … | coef | 94 | -0.3715 | 0.7103 | |
| var_toxic_12_scale ~ area_herb_log_scale + … | coef | 94 | -1.5513 | 0.1208 | |
| select_scale ~ area_herb_log_scale + … | coef | 94 | -0.0253 | 0.9798 | |
| ava_qual_logit_scale ~ area_herb_log_scale + … | coef | 94 | 1.0211 | 0.3072 | |
| on_toxic_logit_scale ~ area_herb_log_scale + … | coef | 94 | 1.1450 | 0.2522 | |
| mean_toxic_conc_scale ~ area_herb_log_scale + … | coef | 94 | 0.5116 | 0.6089 | |
| mean_toxic_conc_scale ~ sl_mean_obs_log_scale + … | coef | 94 | -0.4908 | 0.6236 | |
| mean_toxic_conc_scale ~ var_toxic_12_scale + … | coef | 94 | -1.4198 | 0.1557 | |
| mean_toxic_conc_scale ~ select_scale + … | coef | 94 | 1.1730 | 0.2408 | |
| mean_toxic_conc_scale ~ ava_qual_logit_scale + … | coef | 94 | 0.6698 | 0.5030 |
3.3 Model summary
3.3.1 Model stats
sem_summary$ChiSqNULLsem_summary$AIC AIC AICc K n
1 1198.514 1212.125 63 94sem_summary$R2 Response family link method Marginal Conditional
1 RGR_scale gaussian identity none 0.74 0.76
2 var_toxic_12_scale gaussian identity none 0.72 0.73
3 on_toxic_logit_scale gaussian identity none 0.98 NA
4 mean_toxic_conc_scale gaussian identity none 0.60 NA
5 select_scale gaussian identity none 0.29 NA
6 ava_qual_logit_scale gaussian identity none 0.19 0.21
7 area_herb_log_scale gaussian identity none 0.36 NA
8 sl_mean_obs_log_scale gaussian identity none 0.26 NA3.3.2 Standardized coefficients
knitr::kable(
cbind(sem_summary$coefficients,"arrow_size"=(abs(as.numeric(sem_summary$coefficients$Std.Estimate))*10))
)| Response | Predictor | Estimate | Std.Error | DF | Crit.Value | P.Value | Std.Estimate | star | arrow_size |
|---|---|---|---|---|---|---|---|---|---|
| RGR_scale | cat_pre_wt_log_scale | -3.3671 | 0.4051 | 94 | -8.3124 | 0.0000 | -3.4487000 | *** | 34.4870000 |
| RGR_scale | cat_pre_wt_log_scale_sq | -2.7898 | 0.4174 | 94 | -6.6841 | 0.0000 | -2.7190000 | *** | 27.1900000 |
| RGR_scale | mean_toxic_conc_scale | -0.2339 | 0.0689 | 94 | -3.3935 | 0.0007 | -0.2088000 | *** | 2.0880000 |
| RGR_scale | area_herb_log_scale | 0.7734 | 0.0915 | 94 | 8.4511 | 0.0000 | 0.5289000 | *** | 5.2890000 |
| RGR_scale | sl_mean_obs_log_scale | -0.1649 | 0.0537 | 94 | -3.0704 | 0.0021 | -0.1739000 | ** | 1.7390000 |
| RGR_scale | var_toxic_12_scale | -0.0055 | 0.055 | 94 | -0.0999 | 0.9204 | -0.0061000 | 0.0610000 | |
| RGR_scale | cat_pre_wt_log_scale:area_herb_log_scale | -0.2007 | 0.0912 | 94 | -2.2008 | 0.0277 | -0.1353000 | * | 1.3530000 |
| var_toxic_12_scale | beta_numeric_scale | -0.1088 | 0.0583 | 94 | -1.8663 | 0.0620 | -0.1095000 | 1.0950000 | |
| var_toxic_12_scale | var_high | 1.6971 | 0.115 | 94 | 14.7574 | 0.0000 | 1.7241172 | *** | 17.2411719 |
| var_toxic_12_scale | sl_mean_obs_log_scale | 0.1314 | 0.0647 | 94 | 2.0303 | 0.0423 | 0.1248000 | * | 1.2480000 |
| var_toxic_12_scale | select_scale | -0.2063 | 0.0816 | 94 | -2.5290 | 0.0114 | -0.1643000 | * | 1.6430000 |
| var_toxic_12_scale | ava_qual_logit_scale | -0.2832 | 0.0788 | 94 | -3.5942 | 0.0003 | -0.2702000 | *** | 2.7020000 |
| on_toxic_logit_scale | select_scale | -0.1577 | 0.0194 | 94 | -8.1434 | 0.0000 | -0.1333000 | *** | 1.3330000 |
| on_toxic_logit_scale | ava_qual_logit_scale | -0.9098 | 0.0162 | 94 | -56.2663 | 0.0000 | -0.9207000 | *** | 9.2070000 |
| mean_toxic_conc_scale | beta_numeric_scale | -0.0873 | 0.0533 | 94 | -1.6370 | 0.1016 | -0.1093000 | 1.0930000 | |
| mean_toxic_conc_scale | var_high | -0.8126 | 0.112 | 94 | -7.2565 | 0.0000 | -1.0180231 | *** | 10.1802310 |
| mean_toxic_conc_scale | cat_pre_wt_log_scale | 0.9656 | 0.3702 | 94 | 2.6083 | 0.0091 | 1.1079000 | ** | 11.0790000 |
| mean_toxic_conc_scale | on_toxic_logit_scale | 0.2337 | 0.0609 | 94 | 3.8346 | 0.0001 | 0.2741000 | *** | 2.7410000 |
| mean_toxic_conc_scale | cat_pre_wt_log_scale_sq | 1.0601 | 0.3786 | 94 | 2.7998 | 0.0051 | 1.1574000 | ** | 11.5740000 |
| mean_toxic_conc_scale | var_high:cat_pre_wt_log_scale | -0.2656 | 0.1165 | 94 | -2.2802 | 0.0226 | -0.2404000 | * | 2.4040000 |
| select_scale | beta_numeric_scale | -0.0804 | 0.071 | 94 | -1.1318 | 0.2577 | -0.1016000 | 1.0160000 | |
| select_scale | var_high | -0.0137 | 0.1465 | 94 | -0.0938 | 0.9253 | -0.0173199 | 0.1731994 | |
| select_scale | cat_pre_wt_log_scale | -0.1295 | 0.1203 | 94 | -1.0768 | 0.2816 | -0.1500000 | 1.5000000 | |
| select_scale | ava_qual_logit_scale | 0.3762 | 0.0794 | 94 | 4.7364 | 0.0000 | 0.4506000 | *** | 4.5060000 |
| select_scale | var_high:cat_pre_wt_log_scale | 0.3718 | 0.1538 | 94 | 2.4174 | 0.0156 | 0.3396000 | * | 3.3960000 |
| ava_qual_logit_scale | var_high | 0.4543 | 0.1731 | 94 | 2.6246 | 0.0087 | 0.5111900 | ** | 5.1119004 |
| ava_qual_logit_scale | beta_numeric_scale | 0.1570 | 0.0902 | 94 | 1.7407 | 0.0817 | 0.1657000 | 1.6570000 | |
| ava_qual_logit_scale | cat_pre_wt_log_scale | 0.1843 | 0.1215 | 94 | 1.5176 | 0.1291 | 0.1782000 | 1.7820000 | |
| ava_qual_logit_scale | beta_numeric_scale:cat_pre_wt_log_scale | 0.2136 | 0.0983 | 94 | 2.1721 | 0.0298 | 0.2069000 | * | 2.0690000 |
| area_herb_log_scale | beta_numeric_scale | 0.0016 | 0.0509 | 94 | 0.0319 | 0.9746 | 0.0027000 | 0.0270000 | |
| area_herb_log_scale | var_high | 0.2834 | 0.1046 | 94 | 2.7091 | 0.0067 | 0.4634544 | ** | 4.6345440 |
| area_herb_log_scale | cat_pre_wt_log_scale | 0.5505 | 0.0878 | 94 | 6.2661 | 0.0000 | 0.8245000 | *** | 8.2450000 |
| area_herb_log_scale | var_high:cat_pre_wt_log_scale | -0.3558 | 0.1132 | 94 | -3.1425 | 0.0017 | -0.4204000 | ** | 4.2040000 |
| sl_mean_obs_log_scale | beta_numeric_scale | 0.1876 | 0.0869 | 94 | 2.1586 | 0.0309 | 0.1989000 | * | 1.9890000 |
| sl_mean_obs_log_scale | on_toxic_logit_scale | 0.3847 | 0.1161 | 94 | 3.3120 | 0.0009 | 0.3819000 | *** | 3.8190000 |
| sl_mean_obs_log_scale | var_high | -0.2428 | 0.1768 | 94 | -1.3735 | 0.1696 | -0.2574872 | 2.5748724 | |
| sl_mean_obs_log_scale | select_scale | 0.3570 | 0.1307 | 94 | 2.7319 | 0.0063 | 0.2994000 | ** | 2.9940000 |
| sl_mean_obs_log_scale | cat_pre_wt_log_scale | 1.5178 | 0.5808 | 94 | 2.6131 | 0.0090 | 1.4742000 | ** | 14.7420000 |
| sl_mean_obs_log_scale | cat_pre_wt_log_scale_sq | 1.8475 | 0.612 | 94 | 3.0189 | 0.0025 | 1.7075000 | ** | 17.0750000 |
| ~~var_toxic_12_scale | ~~on_toxic_logit_scale | 0.0823 | - | 94 | 0.7882 | 0.2163 | 0.0823000 | 0.8230000 |
3.3.3 Post hoc contrast
contrast_by_pre_wt(subm_toxin_ingested, "var_high", "trt")$emmeans
cat_pre_wt_log_scale = -2:
var_high emmean SE df lower.CL upper.CL
0 -1.941 0.823 85 -3.577 -0.306
1 -2.223 0.819 85 -3.851 -0.594
cat_pre_wt_log_scale = 2:
var_high emmean SE df lower.CL upper.CL
0 1.921 0.666 85 0.597 3.246
1 0.578 0.629 85 -0.672 1.827
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 -0.281 0.279 85 -1.008 0.3163
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 -1.344 0.236 85 -5.695 <.00011/sd(subm_toxin_ingested$frame$mean_toxic_conc_scale)[1] 1.252797contrast_by_pre_wt(subm_select, "var_high", "trt")$emmeans
cat_pre_wt_log_scale = -2:
var_high emmean SE df lower.CL upper.CL
0 0.251 0.287 86 -0.3197 0.8211
1 -0.507 0.241 86 -0.9854 -0.0278
cat_pre_wt_log_scale = 2:
var_high emmean SE df lower.CL upper.CL
0 -0.267 0.235 86 -0.7336 0.1988
1 0.462 0.201 86 0.0638 0.8610
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 -0.757 0.371 86 -2.041 0.0443
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 0.730 0.307 86 2.375 0.01981/sd(subm_select$frame$select_scale)[1] 1.264229contrast_by_pre_wt(subm_herb, "var_high", "trt")$emmeans
cat_pre_wt_log_scale = -2:
var_high emmean SE df lower.CL upper.CL
0 -1.054 0.208 87 -1.467 -0.641
1 -0.059 0.176 87 -0.410 0.292
cat_pre_wt_log_scale = 2:
var_high emmean SE df lower.CL upper.CL
0 1.148 0.173 87 0.805 1.491
1 0.720 0.143 87 0.436 1.004
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 0.995 0.272 87 3.652 0.0004
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 -0.428 0.224 87 -1.911 0.05931/sd(subm_herb$frame$area_herb_log_scale)[1] 1.635337emmeans::emtrends(subm_rgr,
~ area_herb_log_scale | cat_pre_wt_log_scale,
var = "area_herb_log_scale",
at = list(cat_pre_wt_log_scale = c(-2,2)))cat_pre_wt_log_scale = -2:
area_herb_log_scale area_herb_log_scale.trend SE df lower.CL upper.CL
0.276 1.175 0.190 84 0.7974 1.552
cat_pre_wt_log_scale = 2:
area_herb_log_scale area_herb_log_scale.trend SE df lower.CL upper.CL
0.276 0.372 0.217 84 -0.0603 0.804
Confidence level used: 0.95 sd(subm_rgr$frame$area_herb_log_scale) / sd(subm_rgr$frame$RGR_scale)[1] 0.6838414emmeans::emtrends(subm_ava,
~ beta_numeric_scale | cat_pre_wt_log_scale,
var = "beta_numeric_scale",
at = list(cat_pre_wt_log_scale = c(-2,2)))cat_pre_wt_log_scale = -2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.95e-17 -0.270 0.238 87 -0.744 0.203
cat_pre_wt_log_scale = 2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.95e-17 0.584 0.192 87 0.202 0.966
Results are averaged over the levels of: var_high
Confidence level used: 0.95 sd(subm_ava$frame$beta_numeric_scale) * sd(subm_ava$frame$ava_qual_logit_scale)[1] 0.94757213.4 SEM plots
knitr::include_graphics(paste0(getwd(), "/graphs/manuscript1_figures/SEM_joined.png"), rel_path = FALSE)
3.5 SEM forest plot
knitr::include_graphics(paste0(getwd(),"/graphs/manuscript1_figures/SEM_forest.png"), rel_path = FALSE)
3.6 Misc
3.6.1 Time to pupation
subm_pupt <- glmmTMB(
time_to_pupation ~
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
mean_toxic_conc_scale +
area_herb_log_scale +
sl_mean_obs_log_scale +
var_toxic_12_scale +
(1|session_id),
family = Gamma(link = "log"),
data = d2 %>%
filter(!is.na(time_to_pupation))
); summary(subm_pupt) Family: Gamma ( log )
Formula:
time_to_pupation ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +
mean_toxic_conc_scale + area_herb_log_scale + sl_mean_obs_log_scale +
var_toxic_12_scale + (1 | session_id)
Data: d2 %>% filter(!is.na(time_to_pupation))
AIC BIC logLik deviance df.resid
280.8 302.8 -131.4 262.8 76
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.69e-05 0.006075
Number of obs: 85, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0224
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.205884 0.021976 100.38 < 2e-16 ***
cat_pre_wt_log_scale -1.107727 0.137192 -8.07 6.79e-16 ***
cat_pre_wt_log_scale_sq -0.781595 0.150144 -5.21 1.93e-07 ***
mean_toxic_conc_scale 0.082351 0.024635 3.34 0.000829 ***
area_herb_log_scale -0.139853 0.044549 -3.14 0.001693 **
sl_mean_obs_log_scale 0.058968 0.019041 3.10 0.001956 **
var_toxic_12_scale 0.003529 0.020361 0.17 0.862383
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2(subm_pupt)# R2 for Mixed Models
Conditional R2: 0.888
Marginal R2: 0.8884 Session Info
sessionInfo()R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spat1f_0.0.1 herbivar_0.2.1 forcats_0.5.1 stringr_1.5.0
[5] dplyr_1.1.2 purrr_1.0.1 readr_2.1.2 tidyr_1.3.0
[9] tibble_3.2.1 tidyverse_1.3.1 imager_0.42.13 magrittr_2.0.3
[13] doSNOW_1.0.20 snow_0.4-4 iterators_1.0.14 foreach_1.5.2
[17] glmmTMB_1.1.7 tidybayes_3.0.2 performance_0.10.4 sjPlot_2.8.15
[21] ggpubr_0.6.0 ggplot2_3.4.2 piecewiseSEM_2.3.0 survival_3.3-1
loaded via a namespace (and not attached):
[1] svUnit_1.0.6 splines_4.3.1 later_1.3.0
[4] cellranger_1.1.0 polyclip_1.10-0 reprex_2.0.1
[7] lifecycle_1.0.3 sf_1.0-8 Rdpack_2.3.1
[10] rstatix_0.7.2 doParallel_1.0.17 rprojroot_2.0.3
[13] vroom_1.5.7 processx_3.6.1 lattice_0.21-8
[16] MASS_7.3-60 insight_0.19.3 readbitmap_0.1.5
[19] ggdist_3.1.1 backports_1.4.1 sass_0.4.1
[22] rmarkdown_2.24 jquerylib_0.1.4 yaml_2.3.5
[25] remotes_2.4.2 httpuv_1.6.5 gap_1.2.3-6
[28] sessioninfo_1.2.2 pkgbuild_1.3.1 cowplot_1.1.1
[31] DBI_1.1.3 minqa_1.2.4 RColorBrewer_1.1-3
[34] lubridate_1.8.0 DHARMa_0.4.6 multcomp_1.4-25
[37] abind_1.4-5 pkgload_1.3.2 rvest_1.0.2
[40] imagerExtra_2.0.0 TH.data_1.1-1 tensorA_0.36.2
[43] sandwich_3.0-2 spatstat.utils_3.0-1 units_0.8-0
[46] codetools_0.2-19 MuMIn_1.47.5 xml2_1.3.3
[49] tidyselect_1.2.0 ggeffects_1.2.3 farver_2.1.0
[52] lme4_1.1-29 matrixStats_1.2.0 stats4_4.3.1
[55] jsonlite_1.8.8 e1071_1.7-11 ellipsis_0.3.2
[58] bmp_0.3 emmeans_1.7.5 tools_4.3.1
[61] Rcpp_1.0.8.3 glue_1.6.2 mgcv_1.8-40
[64] xfun_0.39 distributional_0.3.0 usethis_2.1.6
[67] withr_2.5.0 numDeriv_2016.8-1.1 fastmap_1.1.0
[70] boot_1.3-28.1 fansi_1.0.3 callr_3.7.0
[73] digest_0.6.29 R6_2.5.1 mime_0.12
[76] estimability_1.3 colorspace_2.0-3 spatstat.data_3.0-0
[79] jpeg_0.1-9 see_0.7.1 DiagrammeR_1.0.9
[82] utf8_1.2.2 generics_0.1.2 class_7.3-22
[85] prettyunits_1.1.1 httr_1.4.3 htmlwidgets_1.5.4
[88] pkgconfig_2.0.3 gtable_0.3.0 htmltools_0.5.2
[91] carData_3.0-5 profvis_0.3.7 fftwtools_0.9-11
[94] TMB_1.9.4 scales_1.2.1 png_0.1-7
[97] posterior_1.2.2 knitr_1.43 rstudioapi_0.13
[100] tzdb_0.3.0 coda_0.19-4 visNetwork_2.1.0
[103] checkmate_2.1.0 nlme_3.1-162 nloptr_2.0.3
[106] proxy_0.4-27 cachem_1.0.6 zoo_1.8-12
[109] sjlabelled_1.2.0 KernSmooth_2.23-20 parallel_4.3.1
[112] miniUI_0.1.1.1 desc_1.4.1 pillar_1.9.0
[115] grid_4.3.1 vctrs_0.6.3 urlchecker_1.0.1
[118] promises_1.2.0.1 car_3.1-2 dbplyr_2.2.0
[121] arrayhelpers_1.1-0 xtable_1.8-4 evaluate_0.15
[124] amt_0.2.1.0 mvtnorm_1.1-3 cli_3.6.1
[127] compiler_4.3.1 rlang_1.1.1 crayon_1.5.1
[130] ggsignif_0.6.3 labeling_0.4.2 modelr_0.1.8
[133] classInt_0.4-7 RcppArmadillo_0.11.2.0.0 ps_1.7.1
[136] plyr_1.8.7 sjmisc_2.8.9 fs_1.5.2
[139] gap.datasets_0.0.5 stringi_1.7.12 deldir_1.0-6
[142] assertthat_0.2.1 munsell_0.5.0 tiff_0.1-11
[145] spatstat.geom_3.0-3 devtools_2.4.5 bayestestR_0.13.1
[148] Matrix_1.5-4.1 sjstats_0.18.2 hms_1.1.1
[151] bit64_4.0.5 qgam_1.3.4 shiny_1.7.1
[154] highr_0.9 haven_2.5.0 rbibutils_2.2.8
[157] igraph_1.3.2 broom_1.0.5 memoise_2.0.1
[160] bslib_0.3.1 bit_4.0.4 readxl_1.4.0